function Features=getTFNJND(image)
delta=48;
Pa=[0,0,0,0];
Pb=0;
Pc=[0,0,0,0];
[m,n]=size(image);
Matrix=zeros(m,n);
N=zeros(1,42);
deltaP=[-1,-1;-1,0;-1,1;0,1;1,1;1,0;1,-1;0,-1];
for i=1:m
    for j=1:n        
        bg=0;
        for p=1:8
            x=i+deltaP(p,1);
            y=j+deltaP(p,2);
            if(x>m||x<=0||y>n||y<=0)
                Vp=image(i,j);
            else
                Vp=image(x,y);
            end
            bg=bg+Vp;
        end
        bg=(bg+image(i,j))/9;
        JND=0;
        if(bg<=2047)
            JND=17*(1-sqrt(bg/2047))+3;
        else
            JND=3/2048*(bg-2047)+3;
        end
        JND=JND*16;
        delta=0.15*JND;
        
        
        
        Pc=image(i,j);
        if(j-1<=0)
            Pa(1)=Pc;
        else
            Pa(1)=image(i,j-1);
        end
        if(j+1>n)
            Pb(1)=Pc;
        else
            Pb(1)=image(i,j+1);
        end
        if(j+1>n||i-1<=0)
            Pa(3)=Pc;
        else
            Pa(3)=image(i-1,j+1);
        end
        if(j-1<=0||i+1>m)
            Pb(3)=Pc;
        else
            Pb(3)=image(i+1,j-1);
        end
        d1=getTFNnum(Pa(1),Pb(1),Pc,delta);
        d3=getTFNnum(Pa(3),Pb(3),Pc,delta);
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        if(i-1<=0)
            Pa(2)=Pc;
        else
            Pa(2)=image(i-1,j);
        end
        if(i+1>m)
            Pb(2)=Pc;
        else
            Pb(2)=image(i+1,j);
        end
        if(j-1<=0||i-1<=0)
            Pa(4)=Pc;
        else
            Pa(4)=image(i-1,j-1);
        end
        if(j+1>n||i+1>m)
            Pb(4)=Pc;
        else
            Pb(4)=image(i+1,j+1);
        end
        d2=getTFNnum(Pa(2),Pb(2),Pc,delta);
        d4=getTFNnum(Pa(4),Pb(4),Pc,delta);
        tmp=calcTFNAB(d1,d2,d3,d4);
        Matrix(i,j)=tmp;
        N(tmp)=N(tmp)+1;
    end
end
P=N./(m*n);
%Coarseness
Cosrseness=P(42);
%Homogeneity
Homogeneity=P(1);
%Mean Convergence
MeanP=mean(P);
StdP=std(P);
MC=0;
for i=1:42
    MC=MC+abs(i*P(i)-MeanP);
end
MC=MC/StdP;
%Variance
Variance=0;
for i=1:42
    Variance=Variance+P(i)*((i-MeanP)^2);
end
%Entropy
deltaX=2;
deltaY=2;
WTFN=zeros(42,42);
for i=1:m
    for j=1:n
        X=i+deltaX;
        Y=j+deltaY;
        if(X<=m&&X>0&&Y<=n&&Y>0)
            a=Matrix(i,j);
            b=Matrix(X,Y);
            WTFN(a,b)=WTFN(a,b)+1;
        end
    end
end
sumW=sum(sum(WTFN));
PTFN=WTFN./sumW;
EntropyP=entropy(PTFN);
%Similarity
Similarity=sum(sum(PTFN.*PTFN));
%Regularity
PP=zeros(42,42);
for i=1:42
    for j=1:42
        PP(i,j)=1+(i-j)^2;
    end
end
Regularity=sum(sum(PTFN./PP));

Features=[Cosrseness Homogeneity MC Variance EntropyP Similarity Regularity];


